COVID-19 social distancing model

Source: Agents.jl model zoo

using Agents
using Random

Update message: Agents v5
Welcome to this new major version of Agents.jl!
Noteworthy changes:

* Schedulers have been reworked to be more performant and allocate
  less. This means that most scheduler names have been deprecated from
  functions to types, such as `by_type` -> `ByType`. See changelog for
  full list!
* See the CHANGELOG.md or online docs for more!

Let us first create a simple model where balls move around in a continuous space. We need to create agents that comply with ContinuousSpace, i.e. they have a pos and vel fields, both of which are tuples of float numbers.

mutable struct Agent <: AbstractAgent
    id::Int                 # Mandatory Agent identifier
    pos::NTuple{2,Float64}  # Position, required for agents in the ContinuousSpace
    vel::NTuple{2,Float64}  # Moving speeds
    mass::Float64           # Can move or not
end

function ball_model(; speed = 0.002)
    space2d = ContinuousSpace((1, 1), 0.02)
    model = ABM(Agent, space2d, properties = Dict(:dt => 1.0), rng = MersenneTwister(42))

    # Add agents to the model
    for ind in 1:500
        pos = Tuple(rand(model.rng, 2))
        vel = sincos(2π * rand(model.rng)) .* speed
        mass = 1.0
        add_agent!(pos, model, vel, mass)
    end
    return model
end

model = ball_model()

AgentBasedModel with 500 agents of type Agent
 space: periodic continuous space with 50×50 divisions
 scheduler: fastest
 properties: dt

# Agents.move_agent!()
agent_step!(agent, model) = move_agent!(agent, model, model.dt)

agent_step! (generic function with 1 method)

using InteractiveDynamics
using CairoMakie

abmvideo(
    "socialdist1.mp4",
    model,
    agent_step!;
    title = "Ball Model",
    frames = 50,
    spf = 2,
    framerate = 25,
)

using Base64

function display_mp4(filename)
    display("text/html", string("""<video autoplay controls><source src="data:video/x-m4v;base64,""",
        Base64.base64encode(open(read, filename)),"""" type="video/mp4"></video>"""))
end

display_mp4("socialdist1.mp4")

As you can see the agents move in a straight line in a periodic space without interactions. Let’s change that.

Billiard-like interaction

Using the continuous space API:

• interacting_pairs()

• elastic_collision!()

And we redefine the stepping function:

function model_step!(model)
    for (a1, a2) in interacting_pairs(model, 0.012, :nearest)
        elastic_collision!(a1, a2, :mass)
    end
end

model_step! (generic function with 1 method)

model2 = ball_model()

abmvideo(
    "socialdist2.mp4",
    model2,
    agent_step!,
    model_step!;
    title = "Billiard-like",
    frames = 50,
    spf = 2,
    framerate = 25,
)

display_mp4("socialdist2.mp4")

Immovable agents

For the following social distancing example, it will become crucial that some agents don’t move, and can’t be moved (i.e. they stay “isolated”). This is very easy to do with the elastic_collision! function, we only have to make some agents have infinite mass.

model3 = ball_model()

for id in 1:400
    agent = model3[id]
    agent.mass = Inf
    agent.vel = (0.0, 0.0)
end

abmvideo(
    "socialdist3.mp4",
    model3,
    agent_step!,
    model_step!;
    title = "Billiard-like with stationary agents",
    frames = 50,
    spf = 2,
    framerate = 25,
)

display_mp4("socialdist3.mp4")

Adding Virus spread (SIR model)

The agents can be infected with a disease and transfer the disease to other agents around them.

mutable struct PoorSoul <: AbstractAgent
    id::Int
    pos::NTuple{2,Float64}
    vel::NTuple{2,Float64}
    mass::Float64
    days_infected::Int  # number of days since is infected
    status::Symbol  # :S, :I or :R
    β::Float64
end

β is the transmission probability, which we choose to make an agent parameter instead of a model parameter. It reflects the level of hygiene of an individual.

And the model creation function becomes:

const steps_per_day = 24 # One tick per hour

function sir_initiation(;
    infection_period = 30 * steps_per_day,
    detection_time = 14 * steps_per_day,
    reinfection_probability = 0.05,
    isolated = 0.0, # in percentage
    interaction_radius = 0.012,
    dt = 1.0,
    speed = 0.002,
    death_rate = 0.044,
    N = 1000,
    initial_infected = 5,
    seed = 42,
    βmin = 0.4,
    βmax = 0.8,
)

    properties = (;
        infection_period,
        reinfection_probability,
        detection_time,
        death_rate,
        interaction_radius,
        dt,
    )
    space = ContinuousSpace((1,1), 0.02)
    model = ABM(PoorSoul, space, properties = Dict(pairs(properties)), rng = MersenneTwister(seed))

    # Add initial individual agents
    for ind in 1:N
        pos = Tuple(rand(model.rng, 2))
        status = ind ≤ N - initial_infected ? :S : :I
        isisolated = ind ≤ isolated * N
        mass = isisolated ? Inf : 1.0
        vel = isisolated ? (0.0, 0.0) : sincos(2π * rand(model.rng)) .* speed

        β = (βmax - βmin) * rand(model.rng) + βmin
        add_agent!(pos, model, vel, mass, 0, status, β)
    end

    return model
end

sir_initiation (generic function with 1 method)

To visualize this model, we will use black color for the susceptible, red for the infected infected and green for the recovered.

sir_model = sir_initiation()

sir_colors(a) = a.status == :S ? "#2b2b33" : a.status == :I ? "#bf2642" : "#338c54"

fig, abmstepper = abmplot(sir_model; ac = sir_colors)
fig # display figure

Modify the model_step! function to simulate transmission.

function transmit!(a1, a2, rp)

    # for transmission, only 1 can have the disease (otherwise nothing happens)
    count(a.status == :I for a in (a1, a2)) ≠ 1 && return
    infected, healthy = a1.status == :I ? (a1, a2) : (a2, a1)

    # Lucky and not infected
    rand(model.rng) > infected.β && return

    # Risk of reinfection
   if healthy.status == :R
        rand(model.rng) > rp && return
    end

    # You got virus
    healthy.status = :I
end

function sir_model_step!(model)
    r = model.interaction_radius
    for (a1, a2) in interacting_pairs(model, r, :all)
        transmit!(a1, a2, model.reinfection_probability)
        elastic_collision!(a1, a2, :mass)
    end
end

# Agent-specific functions
function update!(agent) 
    if agent.status == :I
        agent.days_infected += 1
    end
end

function recover_or_die!(agent, model)
    if agent.days_infected ≥ model.infection_period
        if rand(model.rng) ≤ model.death_rate
            kill_agent!(agent, model)
        else
            agent.status = :R
            agent.days_infected = 0
        end
    end
end

function sir_agent_step!(agent, model)
    move_agent!(agent, model, model.dt)
    update!(agent)
    recover_or_die!(agent, model)
end

sir_agent_step! (generic function with 1 method)

Running with default parameters.

sir_model = sir_initiation()

abmvideo(
    "socialdist4.mp4",
    sir_model,
    sir_agent_step!,
    sir_model_step!;
    title = "SIR model",
    frames = 200,
    ac = sir_colors,
    as = 10,
    spf = 2,
    framerate = 20,
)

display_mp4("socialdist4.mp4")

Analyzing exponential spread

infected(x) = count(i == :I for i in x)
recovered(x) = count(i == :R for i in x)
# Aggregated data for number of infected and recovered indivisuals
adata = [(:status, infected), (:status, recovered)]

2-element Vector{Tuple{Symbol, Function}}:
 (:status, infected)
 (:status, recovered)

# Try different parameters
r1, r2 = 0.02, 0.05
β1, β2 = 0.5, 0.1
sir_model1 = sir_initiation(reinfection_probability = r1, βmax = β1)
sir_model2 = sir_initiation(reinfection_probability = r2, βmax = β1)
sir_model3 = sir_initiation(reinfection_probability = r1, βmax = β2)

data1, _ = run!(sir_model1, sir_agent_step!, sir_model_step!, 3000; adata)
data2, _ = run!(sir_model2, sir_agent_step!, sir_model_step!, 3000; adata)
data3, _ = run!(sir_model3, sir_agent_step!, sir_model_step!, 3000; adata)

data1[(end-10):end, :]

11 rows × 3 columns

	step	infected_status	recovered_status
	Int64	Int64	Int64
1	2990	76	875
2	2991	76	875
3	2992	76	875
4	2993	77	874
5	2994	77	874
6	2995	77	874
7	2996	77	874
8	2997	77	874
9	2998	77	874
10	2999	77	874
11	3000	77	874

using CairoMakie

figure = Figure()
ax = figure[1, 1] = Axis(figure; ylabel = "Infected", xlabel="Steps")
l1 = lines!(ax, data1[:, dataname((:status, infected))], color = :orange)
l2 = lines!(ax, data2[:, dataname((:status, infected))], color = :blue)
l3 = lines!(ax, data3[:, dataname((:status, infected))], color = :green)
figure[1, 2] = Legend(figure, [l1, l2, l3], ["r=$r1, beta=$β1", "r=$r2, beta=$β1", "r=$r1, beta=$β2"])

figure

Social distancing

The best way to model social distancing is to make some agents simply not move (which feels like it approximates reality better).

sir_model = sir_initiation(isolated = 0.85)
abmvideo(
    "socialdist5.mp4",
    sir_model,
    sir_agent_step!,
    sir_model_step!;
    title = "Social Distancing",
    frames = 200,
    spf = 2,
    ac = sir_colors,
    framerate = 20,
)

display_mp4("socialdist5.mp4")

r4 = 0.02
sir_model4 = sir_initiation(reinfection_probability = r4, βmax = β1, isolated = 0.85)

data4, _ = run!(sir_model4, sir_agent_step!, sir_model_step!, 3000; adata)

l4 = lines!(ax, data4[:, dataname((:status, infected))], color = :red)
figure[1, 2] = Legend(
    figure,
    [l1, l2, l3, l4],
    ["r=$r1, beta=$β1", "r=$r2, beta=$β1", "r=$r1, beta=$β2", "r=$r4, social distancing"],
)
figure